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ABSTRACT 


Accurate identification of fish and fish products, from eggs to adults, is important in many areas. Grey mullets 
of the family Mugilidae are distributed worldwide and inhabit marine, estuarine, and freshwater environments 
in all tropical and temperate regions. Various Mugilid species are commercially important species in fishery 
and aquaculture of many countries. For the present study we have chosen two Mugilid genes with different 
phylogenetic signals: relatively variable mitochondrial cytochrome oxidase subunit I (COI) and conservative 
nuclear rhodopsin (RHO). We examined their diversity within and among 9 Mugilid species belonging to 4 
genera, many of which have been examined from multiple specimens, with the goal of determining whether 
DNA barcoding can achieve unambiguous species recognition of Mugilid species. The data obtained showed 
that information based on COI sequences was diagnostic not only for species-level identification but also for 
recognition of intraspecific units, e.g., allopatric populations of circumtropical Mugil cephalus, or even native 
and acclimatized specimens of Chelon haematocheila. All RHO sequences appeared strictly species specific. 
Based on the data obtained, we conclude that COI, as well as RHO sequencing can be used to unambiguously 
identify fish species. Topologies of phylogeny based on RHO and COI sequences coincided with each other, 
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while together they had a good phylogenetic signal. 
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INTRODUCTION 


The unequivocal identification and classification of living 
organisms to the species level frequently relies on genetic 
evidence. Specific DNA sequences act as unrepeatable signa- 
tures and therefore constitute a unique DNA barcode for each 
species. Hebert et al. (2003) proposed that a single gene 
sequence would be sufficient to differentiate all, or at least 
the vast majority of animal species, and proposed the use of 
the mitochondrial DNA gene cytochrome oxidase subunit I 
(COD) as a global bioidentification system for animals. Ini- 
tiatives, such as The Barcode of Life Database (http: //www. 
barcodinglife.org) including The Fish Barcode of Life (http: 
//www.fishbolL.org), use a DNA-based identification system 
based on a relatively small fragment of COI. 

However, the approach is not without controversy (Lips- 


comb et al., 2003; Moritz and Cicero, 2004). For a barcoding 
approach to species identification to succeed, within-species 
DNA sequences need to be more similar to one another than 
to sequences in different species. Recent studies show that 
this is generally the case, but there are exceptions. Hybridiza- 
tion among species would create taxonomic uncertainty: mito- 
chondrial DNA is maternally inherited and any hybrid or 
subsequent generation would have the maternal species DNA 
only. Thus, the idea of a multi-locus DNA barcoding appro- 
ach is progressively emerging, and limitations of mtDNA 
underline the requirement of nuclear regions (Frézal and Le- 
blois, 2008). It is clear that longer length DNA barcodes will 
provide more efficient identification labels. Barcode effici- 
ency can be further improved by the simultaneous use of 
genes showing different evolutionary rates and genomic posi- 
tions. 
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RHOopsin, a nuclear gene encoding a transmembrane G- 
protein-coupled receptor for visual transduction cascade, pro- 
ved to be very suitable for barcoding purposes because it has 
a single copy in the genome, does not contain introns in Tel- 
eostei fishes, its evolutionary rate in fish is less than 2-fold 
lower than, for example, cyt b gene (0.167 vs. 0.247), and 
RHO trees contain a high number of well supported Teleostei 
clades (Chen et al., 2003). 

Here we examine the diversity of two genes with different 
phylogenetic signal: relatively variable mitochondrial COI 
and conservative nuclear RHOopsin (RHO) within and among 
9 Mugilid species belonging to 4 genera, many of which have 
been examined from multiple specimens, with the goal of 
determining whether DNA barcoding can achieve unambigu- 
ous species recognition of Mugilid species. 

Grey mullets of the family Mugilidae (Pisces, Mugilifor- 
mes) are distributed worldwide and inhabit marine, estuarine, 
and freshwater environments in all tropical and temperate 
regions. Various Mugilid species are commercially important 
species in fishery and aquaculture of many countries. Mugi- 
lid taxonomy still has not been finalized. The family includes 
from 14 to 20 genera recognized as valid according to the 
most recent revisions and 55 currently recognized species 
(Thomson, 1997; Nelson, 2006). Most of them are represen- 
tatives of Liza, Chelon, and Mugil genera. 


MATERIALS AND METHODS 


Sampling, DNA extraction, polymerase chain reaction 
(PCR) amplification and sequencing 

One to seven specimens were collected for each of the nine 
species of Mugilidae (Table 1). DNA extracts were prepared 
from either muscle or heart tissue or fin clips preserved in 
95% ethanol alcohol according to the protocol of Sambrook 
et al. (1989). Amplification of the COI gene fragments was 
carried out using the primers: FishF1 (5’-TCAACCAACCA 
CAAAGACATTGGCAC-3’) and FishR 1 (5'-TAGACTTCT 
GGGTGGCCAAAGAATCA-3’) (Ward et al., 1994). Primers 
RHO0545 (5'-GCAAGCCCATCAGCAACTTCCG-3’) and 
RHOO1039 (5'-TGCTTGTTCATGCAGATGTAGA-3’) 
(Chen et al., 2003) were used to amplify the RHO gene frag- 
ments. 

PCRs were conducted for both genes with a total volume 
of 25 uL consisting of approximately 50 ng of template DNA, 
0.25 uM of each primer, 2.5 uL of 10 x PCR reaction buffer, 
2mM of each dNTP, and 1 U of Taq DNA polymerase (Sib- 
enzyme, Russia). PCR consisted of an initial denaturation 
step at 94°C for 5 min., followed by 30 cycles of denatura- 
tion at 94°C for 30s, annealing for 30 s at 60°C for COI gene 
fragment and 55°C for RHO gene fragments and extension 
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at 72°C for 1 min. The terminal extension was at 72°C for 10 
min. Aliquots (3 uL) of amplicons were examined in 1.5% 
gels, stained with ethidium bromide, and photographed upon 
transillumination. DNA fragments were cut from the agarose 
gels, frozen-thawed, and re-PCRed under the same condi- 
tions. The rePCRed products were precipitated with ethyl 
alcohol and vacuum dried. Then they were sequenced under 
conditions recommended by the manufacturer, using the 
BigDye Terminator v.3.1 Cycle Sequencing Kit (Applied 
Biosystems, Foster City, CA, USA). Electrophoreses of sequ- 
encing products were made in ABI Prism 3130 DNA sequ- 
encer in 50 cm capillary array with polymer POP-7. 


Phylogenetic analyses 

The sequences of COI and RHO gene fragments were ana- 
lyzed both separately and together. DNA sequences were 
aligned using Clustal W 1.8 in MEGA 5.05. Evolution mod- 
els and parameters were separately selected for COI and 
RHO and were then used in the following analyses. 

COI sequence divergences were calculated using the Hase- 
gawa, Kishino, and Yano (HKY85+ G) model (Hasegawa 
et al., 1985) as the appropriate model of sequence evolution 
chosen on the basis of hierarchical likelihood ratio tests 
(LRTs) as implemented in ModelTest 3.8 (Posada and Cran- 
dall, 1998). RHO sequence divergences were estimated using 
MEGA 5.05 by Kimura 2-parameter model. Saturation levels 
of each gene fragment and the combined data set were asse- 
ssed by plotting transitions and transvertions accumulation 
for each pair of haplotypes against Hasegawa, Kishino, and 
Yano distances for COI and Kimura 2 distances for RHO in 
DAMBE 4.5.2 and MEGA 5.05. 

Data were analyzed using three approaches, neighbour 
joining (NJ), maximum parsimony (MP), and Bayesian infer- 
ence (BI). NJ, MP and BI trees were created to provide a 
graphic representation of the patterning of divergence bet- 
ween species computed in MEGA ver. 5.05 and MrBayes 
3.2.1. MP analyses were performed using heuristic searches 
with 50 random stepwise additions and tree bisection-recon- 
nection branch swapping. Bootstrap analyses were used to 
assess the relative robustness of branches with 1,000 repli- 
cates. BI analyses were conducted in MrBayes 3.2.1 with the 
selected best fit models and parameters. Each BI analysis was 
run over 4,000,000 generations using four Markov Chain 
Monte Carlo chains and every 200th tree saved. Four thou- 
sand trees were discarded as burn-in. The robustness of trees 
was tested using posterior probability. 

Scomberomorus cavalla (Liu, 2010) and Tribolodon hako- 
nensis were included in all analyses as outgroups. The sequ- 
ences of Scomberomorus cavalla were taken from GenBank 
(Table 1). 

All sequences are available in Genbank under the accession 
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Table 1. List of Mugilidae species studied with collection locality, number of individuals sequenced, GenBank accession number 


GenBank accession No. 


Species, F . Individual 
Collection locality 
common name sequenced CO1 rhod 
Chelon labrosus, Tyrrhenian Sea, Orbetello Lagoon, Italy 2 EU392233 JX298796 
thicklip mullet 
Liza aurata, Azov Sea, Tuzla, Cape, Russia 2 EU392234 JX298792 
golden grey mullet Tyrrhenian Sea, Orbetello Lagoon, Italy 2 EU392235 JX298791 
Liza haematocheila, Azov Sea, Tuzla Cape, Russia 2 EU392236 JX298799 
redlip mullet East Sea, South Primorye, Vostok Bay, Russia 2 EU392237 JX298800 
Bioresearch Center, Academia Sinica, L EU392238 JX298801 
ID ASIZP0802104, Taiwan 
Chelon macrolepis, Bioresearch Center, Academia Sinica, 2 EU392239 JX298798 
largescale mullet ID ASIZP0800182, Taiwan 
Liza ramado, Tyrrhenian Sea, Orbetello Lagoon, Italy 2 EU392240 JX298797 
thinlip mullet 
Liza saliens, Tyrrhenian Sea, Orbetello Lagoon, Italy L EU392241 JX298795 
leaping mullet 
Chelon subviridis, Bioresearch Center, Academia Sinica, 1 EU392242 JX298794 
greenback mullet ID ASIZP0800195, Taiwan 
Mugil cephalus, Azov Sea, Tuzla Cape, Russia 2 EU392243 JX298802 
flathead mullet East Sea, South Primorye, Vostok Bay Russia 2 EU392244 JX298803 
Mediterranean Sea, Sardinia, 2 EU392245 JX298804 
San Giovanni Lagoon, Italy 
Bioresearch Center, Academia Sinica, 1 EU392246 JX298805 
ID ASIZP0800468, Taiwan 
Moolgarda cunnesius, Bioresearch Center, Academia Sinica, 1 EU392247 JX298793 
longarm mullet ID ASIZP0800182, Taiwan 
Scomberomorus cavalla GenBank 1 GU225662 DQ874799 
Tribolodon hakonensis East Sea, South Primorye, Vostok Bay, Russia 1 EU392225 GU479895 


numbers shown in the Table 1. 


RESULTS 


Sequence characteristics 

A total of 9 Mugilid species were analyzed, generating (be- 
cause of multiple specimens for most species) a total of 60 
sequences. Approximately 652 bp were amplified and sequ- 
enced from the 5' region of the COI gene from mitochondrial 
DNA. No insertions, deletions or stop codons were observed 
in any sequence. One hundred ninty-eight variable sites were 
revealed, most of them (88.9%) were attributable to the the 
3rd codon base, 10.1% to the 2nd, and 1% to the Ist codon 
base. Fourteen nucleotide substitutions resulted in amino 
acid replacement. 

The GC content of the Liza-Chelon-Moolgarda species 
was higher of M. cephalus (48.1% vs. 45.0%) due to a higher 
GC content of the first (58% vs. 55%) and, especially, second 
(54% vs. 43%) and third (44% vs. 32%) codon base in the 
former. 
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The HK Y85 distance of individuals within L. haematoch- 
eila was 0.17% compared to 3.22% within M. cephalus speci- 
mens. Mean divergence among congeneric species was 
12.54%. Mean divergence among species within the family 
increases to 21.8% Mugil vs. Liza, and Mugil vs. Chelon, 
20.5% Mugil vs. Moolgarda, 18.2% Moolgarda vs. Liza, and 
Moolgarda vs. Chelon, and only 10.6% Liza vs. Chelon. 

RHOopsin sequences were 515 bp length. All RHOopsin 
sequences were species-specific. We observed no intraspe- 
cific variability, except one nucleotide substitution in 3rd 
codon position in L. aurata individuals from Azov and Medi- 
terranean Seas that did not result in amino acid replacement. 

Forty-four variable sites were revealed, 28 of them were 
found to be parsimony-informative, most of them (67%) were 
attributable to the 3rd codon base, 5% to the 2nd, and 28% 
to the Ist codon base. Eight nucleotide substitutions resulted 
in amino acid replacement. GC content was about the same 
in all the species studied-53-55%. The level of divergence 
due to conservativeness of RHO gene was much lower com- 
pared to COI ranging from only 0.4% between L. saliens and 
C. labrosus and 5% between L. aurata and M. cephalus. 
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Fig. 1. Bayesian inference (BI) tree based on combined cytochrome oxidase subunit I and RHO gene fragments. Numbers above 
the branches represent the bootstrap values >75 for neighbour joining and maximum parsimony trees as well as posterior 
probabilities values for BI-tree with virtually the same topology. A, Azov Sea; M, Mediterranean Sea; J, East Sea; T, Taiwan. 


The phylogenetic relationships of Mugilidae 

Although DNA barcoding aims to develop species identifica- 
tion systems, some phylogenetic signal was apparent in the 
data. No saturation was detected for each dataset of COI and 
RHO sequences. The plots of transversion and transitions 
against genetic distances were almost linear. Therefore, all 
substitutions in COI and RHO gene fragments were used 
for phylogenetic reconstructions. 

In all cases, no significant differences among the NJ, MP, 
and BI tree topologies based on COI and RHO gene sequ- 
ences were found. The consensus tree obtained from Baye- 
sian analysis of COI dataset was completely identical in 
topology to the one produced using NJ and MP analysis. The 
topologies of NJ, MP, and BI trees based on RHO gene 
analyses were very similar to each other. Thus only the Baye- 
sian tree based on combined COI and RHO gene fragments 
was presented together with major nodal support values for 
NJ and MP (Fig. 1). 

Haplotypes of the same species were always clustered 
together in both phylogenetic reconstructions based on COI, 
RHO and combined COI and RHO gene fragments dataset. 
All haplotypes of M. cephalus sampled from the seas of Azov, 
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Mediterranean and Japan and near the water of Taiwan form- 
ed the monophyletic cluster on the phylogenetic tree based 
on COI dataset as well as RHO sequences with high bootstrap 
and posterior probabilities values (BP, PP). The haplotypes 
of M. cephalus appeared basal as sister group to all other 
Mugilidae on both phylogenetic trees (Fig. 1). L. haemato- 
cheila haplotypes sampled in the East Sea and Okhotsk and 
near Taiwan were clustered together with high nodal support. 
The closest relative of L. haematocheila was C. subviridis 
(Fig. 1). Chelon macrolepis was intermediate between Medi- 
terranean mullets and L. haematocheila-C. subviridis cluster. 
Chelon labrosus grouped with L. aurata, L. ramada, L. sali- 
ens to form a monophyletic clade (Fig. 1), which turned out 
to exclusively comprise species distributed in Mediterranean- 
Azov waters. Liza aurata haplotypes sampled in Mediterran- 
ean and Azov seas clustered together (Fig. 1). The other 
Chelon species C. subviridis together with L. haematochalia 
sampled, all from the Pacific as well as the specimens accli- 
matized on Azov Sea from the East Sea, formed a distinct 
clade. 

Although the phylogenetic trees based on COI and RHO 
sequences were significantly similar to each other, there were 
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differences between them. Haplotypes of M. cephalus sam- 
pled in Azov-Mediterranean, the East Sea and Taiwan resolv- 
ed into three clusters on the phylogenetic tree based on COI 
sequences with high BP and PP support. Haplotypes of L. 
haematocheila consisted of three lineages, one of which was 
acclimatized individuals sampled in Azov Sea, and two were 
from the East Sea and near Taiwan. No variation on RHO 
sequences was revealed within species M. cephalus and L. 
haematoheila. Topology of C. macrolepis and M. cunnesius 
branches was unstable in the phylogenetic tree based on RHO 
sequences, BP and PP support were low. 


DISCUSSION 


Molecular markers are used for explanation and refining the 
taxonomic and phylogenetic relationships in the groups of 
the species. Furthermore, molecular markers are used for the 
accurate and unambiguous identification of fish and fish pro- 
duct from eggs to adult, and it is important in many areas. 

Although phylogenetic analyses of Mugilidae species were 
carried out based on 16S and 12S rRNA mtDNA sequences 
(Liu et al., 2010), Phe 12SrRNA, cyt b and COI (Heras et al., 
2009), three mtDNA loci 16SrRNA, COI and cyt b (Durand 
et al., 2012), the author supposed that employing slower- 
evolved and independent nuclear genes might prove helpful 
to provide molecular systematic of Mugilidae. 

For the barcoding we proposed to use not only variable 
mitochondrial cytochrome oxidase subunit I, but also the 
highly conserved nuclear RHOopsin gene, which has not pre- 
viously been used for such purposes. The use of the nuclear 
genes for the barcoding is especially important in the case of 
hybridization among species. The mitochondrial COI gene 
is the most efficient gene marker. However, in some cases 
5'-COI “barcoding” region is not enough and even be mis- 
leading for purposes of identification due interspecific COI 
recombination reflecting putative historical hybridization 
events between species (Balakirev et al., 2012). 

All 9 species can be differentiated by both COI and RHO. 
The lack of stop codons is consistent with all amplified sequ- 
ences being functional mitochondrial COI sequences, and 
that, together with the fact that all amplified sequences were 
about 628 -652 bp in length, suggests that NUMTs (nuclear 
DNA sequences originating from mitochondrial DNA sequ- 
ences) were not sequenced (vertebrate NUMTS are typically 
smaller than 600 bp (Zhang and Hewitt, 1996). 

The data obtained showed that information based on COI 
sequences was diagnostic not only for species level identifi- 
cation but also for recognition of intraspecific units, e.g., 
allopatric populations of circumtropical M. cephalus, whose 
pronounced population genetic structure is well-known, or 
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even native and acclimatized specimens of L. haematocheila. 

COI sequences of each specimen of each sample proved 
to be identical. We found no difference within each of M. 
cephalus and L. aurata samples from Azov and Mediterran- 
ean Seas. Sequences of geographically remote L. haemato- 
cheila (Azov Sea, East Sea, and Taiwan) and M. cephalus 
(Azov-Mediterranean, East Sea, and Taiwan) had their own 
haplotypes, 1-2 nucleotide substitutions being observed 
between the redlip mullet specimens causing | amino acid 
substitution, while 17-31 nucleotide substitutions being re- 
vealed between the flathead mullet specimens, which result- 
ed in 4 amino acid substitutions. COI short DNA sequence 
provides sufficient identification labels in terms of nucleotide 
positions to discriminate even between congeneric fish spec- 
ies (Hebert et al., 2003) despite the intraspecific variation 
that was especially pronounced between geographically dis- 
tant M. cephalus specimens. 

All RHOopsin sequences were species-specific. Even M. 
cephalus specimens that were geographically remote and 
deeply diverged by mtDNA appeared to be identical. Thus, 
RHO sequences can be considered useful for barcoding pur- 
poses owing to low or zero intraspecific variation and the 
absence of problems connected with introgressive hybridi- 
zation as may emerge in the case of mtDNA. 

The clades revealed after bootstrapping generally corre- 
sponded well with expectations. Topologies of NJ, MP and 
BI trees conformed to our previous data (Semina et al., 2007) 
based on PCR-restiction fragment length polymorphism- 
analysis of extended mtDNA segments. 

The data obtained indicated that M. cephalus was the most 
genetically distant from the rest of the species, which con- 
formed to the results of other authors based on mitochondr- 
ial DNA, allozyme, gemoglobin, and karyologic study (Cat- 
audella et al., 1974; Rizotti, 1993; Caldara et al., 1996; Rossi 
et al., 1998, 2004; Gornung et al., 2001, 2007; Papasotiro- 
poulos et al., 2001, 2002, 2007; Turan et al., 2005; Imsiridou 
et al., 2007). Three mitochondrial lineages of M. cephalus 
based on COI sequences had the regional distribution (Fig. 
1). One corresponded the Mediterranean lineage included 
M. cephalus from Morocco, Tunisia and France and two con- 
formed to the lineages revealed in the Pacific region (Durand 
et al., 2012). We found that all Mediterranean-Azov Chelon 
and Liza species clustered together as well as all Pacific 
Chelon and Liza. The results conformed the close genetic 
relation of Chelon and Liza representatives. Many other re- 
searchers also stated the close relationships of these genera. 
Morphologic (pharyngobranchial organ), cytogenetic, genet- 
ic (mitochondrial DNA), allozyme analyses performed by 
other authors, revealed no substantial difference between 
Mediterranean Liza species and C. labrosus (Cataudella et 
al., 1974; Caldara et al., 1996; Gornung et al., 2001, 2007; 
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Murgia et al., 2002, Papasotiropoulos et al., 2002, 2007; 
Rossi et al., 2004; Turan et al., 2005; Imsiridou et al., 2007; 
Durand et al., 2012). We consider that molecular phylogeny 
offered no evidence for distinguishing Liza from Chelon, 
leading us to challenge their validity. Therefore, taking into 
account the scientific issue data, our results point out the 
necessity for taxonomic revision of Chelon and Liza genera. 

Based on the data obtained, we conclude that COI, as well 
as RHO sequencing, can be used to unambiguously identify 
fish species. Topologies of phylogeny based on COI and 
RHO sequences coincided with each other, while together 
they had a good phylogenetic signal. 
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